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Abstract 

Finite  difference  methods  for  solving  problems  of  time-harmonic 
acoustics  are  developed  and  analyzed.  Multi-dimensional  inhomoge¬ 
neous  problems  with  variable,  possibly  discontinuous,  coefficients  are 
considered,  accounting  for  the  effects  of  employing  non-uniform  grids. 
A  weighted-average  representation  is  less  sensitive  to  transition  in 
wave  resolution  (due  to  variable  wave  numbers  or  non-uniform  grids) 
than  the  standard  pointwise  representation.  Further  enhancement  in 
method  performance  is  obtmned  by  basing  the  stencils  on  generaliza¬ 
tions  of  Fade  approximation,  or  geneialized  definitions  of  the  deriva¬ 
tive,  reducing  spurious  di.spersion,  anisotropy  and  reflection,  and  by 
improving  the  representation  of  source  terms.  The  resulting  schemes 
have  fourth-order  accuiate  local  truncation  error  on  uniform  grids  and 
third  order  in  the  non-uniform  case.  Guidelines  for  discretization  per¬ 
taining  to  grid  orientation  and  resolution  are  presented. 
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1  Introduction 


Boundary-value  problems  governed  by  the  Helmholtz  equation  model  propa¬ 
gation  and  evanescence  of  time-harmonic  waves,  describing  a  variety  of  physi¬ 
cal  phenomena,  including  acoustic,  elastic  and  electromagnetic  waves.  When 
the  wavelength  is  of  the  same  order  as  characteristic  length  scales  cisymptotic 
methods  usually  cannot  be  employed  and  standard  numerical  methods  such 
^ls  boundary  element,  finite  element  or  finite  difference  methods  are  often 
required. 

Finite  difference  methods  are  a  prevalent  comjuilational  technique  that 
applies  to  variable  coefficient  as  well  as  nonlinear  problems.  A  general  frame¬ 
work  for  deriving  higher-order  finite  difference  schemes  Wtis  proposed  by 
Lynch  and  Rice  for  ordinary  differential  equations  fl]  and  elliptic  partial 
differential  equations  [2],  and  applied  to  the  Helmholtz  equation  [3].  The 
coefficients  of  the  stencil  in  this  method  are  computed  by  solving  a  local  sys¬ 
tem  of  equations  so  that  it  is  exact  on  a  given  space  of  polynomials.  This  is 
repeated  at  every  gird  point  at  which  the  solution  is  unknown,  contributing 
to  the  cost  of  computation.  Accuracy  is  further  enhanced  by  judiciously  se¬ 
lecting  the  points  at  which  source  terms  are  evaluated  and  computing  their 
coefficients  in  the  same  way. 

A  family  of  finite  difference  schemes  with  improved  representation  of  a 
range  of  wave  numbers  is  presented  and  analyzed  in  [4].  Tam  and  Webb  [5] 
optimize  the  dispersion  properf.es  of  higher-order  finite  difference  discretiza¬ 
tion  of  the  linearized  Euler  equations.  In  this  approach  the  order  of  accuracy 
of  a  stencil  is  allowed  tc  drop,  freeing  a  parameter  that  is  then  designed  to 
improve  dispersion  respon.s  '. 

Finite  difference  equations  can  be  obtained  by  replacing  the  limit  that 
defines  the  derivative  with  a  finite  counterpart.  The  familiar  definition  of 
the  derivative  may  be  generalized  by  introducing  a  re.solution-dependent  pa¬ 
rameter  leading  to  improved  performance  of  the  tliscrete  methods.  As  long  as 
the  limit  behavior  is  unaltered  consistency  of  the  approximation  is  retained. 
This  idea  was  introduced  by  Mickens  and  employed  a.s  a  means  of  generat¬ 
ing  stable  finite  difference  schemes  on  uniform  grids  for  several  differential 
equations  in  one  spatial  dimension  ((6,  7]  and  references  therein).  Similar 
discrete  equations  are  obtained  by  new  classes  of  finite  element  methods  for 
a  variety  of  applications,  including  wave  piopagation  (e.g.,  [8]  and  references 
therein).  It  should  be  noted  that  analysis  of  the  computation  ot  waves  [9] 
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indicates  that  accuracy  depends  not  only  on  the  product  of  wave  nmnbcr 
and  grid  size  (related  to  resolution),  but  also  on  product  of  other  powers  of 
these  quantities. 

In  this  work  we  apply  generalizations  of  the  definition  of  the  derivative  and 
of  Fade  approximation,  to  finite  difference  stencils  for  the  Helmholtz  equa¬ 
tion  in  order  to  obtain  improved  dispersion  behavior.  Contrary  to  HODIE 
methods,  the  coefficients  are  obtained  explicitly.  Multi-dimensional  inhomo¬ 
geneous  problems  with  variable,  possibly  discontinuous,  coefficients  are  con¬ 
sidered,  accounting  for  the  effects  of  employing  non-uniform  grids.  Several 
finite  difference  stencils  in  one  and  two  dimensions  are  presented  in  Sec.  2. 
The  analysis  of  the  numerical  methods  gradually  builds  up  to  the  general 
case.  Performance  of  the  various  formulations  for  homogeneous  problems 
with  constant  coefficients  on  uniform  grids  is  examined  by  dispersion  anal¬ 
ysis  in  Sec.  3.  This  tool  is  employed  to  define  the  resolution-dependent 
parameter  for  improved  performance.  In  Sec.  4  the  effect  of  the  direction 
of  wave  propagation  relative  to  grid  lines  is  accounted  for.  The  effects  of 
non-uniform  grids  and  discontinuities  in  physical  properties  are  investigated 
in  Sec.  5.  Standard  finite  difference  methods  are  often  not  well-suited  for 
interface  problems  (see,  e.g.,  [10,  pp.  17-21]).  However,  appropriate  repre¬ 
sentation  preserves  the  order  of  accuracy  of  the  continuous-coefficient,  and 
even  constant-coefficient  case.  (Issues  related  to  curved  interfaces,  as  well  as 
curved  domain  boundaries  are  not  treated  herein.)  The  results  of  these  anal¬ 
yses  are  corroborated  by  means  of  local  truncation  error  analysis  in  Sec.  6, 
accounting  also  for  the  effects  of  source  terms.  Numerical  testing  of  these 
stencils  will  be  peiformed  in  future  work. 

2  Discrete  Formulations 

The  Helmholtz  equation  is 


-t-  <(>  -f-  /  =  0  ( 1 ) 

where  k  =  ujfco  is  the  wave  number,  u>  is  the  angular  frequency  and  co  is  the 
speed  of  sound,  and  /  is  a  given  source  distribution.  Although  not  explic  itly 
addressed  in  the  following,  the  case  of  A:*  <  0  which  rorrc'sponds  to  evanescent 
waves  or  singula-  diffusion  problems  is  also  covered.  An  inhomogeneous 
medium  is  represented  by  spatial  variations  in  k(x). 


2.1  One  dimension 


Consider  a  uniform  grid  of  si2e  h  with  points  at  xj  =  jh.  A  typical  start¬ 
ing  point  is  baised  on  the  standard  finite  difference  stencil  for  the  second 
derivative 


Dxr<l>j 


(2) 


and  pointwise  (PT)  representation  of  undifferentiated  terms 


+  k^(i>j  +  /,  =  0 


(3) 


where  is  the  discrete  solution  at  point  j  and  fj  =  fixj).  On  a  non-uniform 
grid  this  generalizes  to 


f  <f>m  -  <i>j 

V 


<()j  -  4>j-i^  I k 


+  fcVi  +  fj  =  o 


(4) 


where  h~  and  h'*'  are  the  grid  size  before  and  after  point  j,  respectively.  For 
a  discontinuity  in  physical  properties  at  point  j  the  stencil  becomes 


<f>}  <t>i 


/i+ 


/f+  +  h- 


^6,+/,  =0  (5) 


where  k~  and  k'^  are  the  wave  numbers  before  and  after  point  ji,  respec¬ 
tively.  These  may  be  also  considered  as  piecewise-constant  approximations 
of  variable  coefficients. 

The  undifferentiated  terms  may  be  represented  by  a  weighted  average 
(wa)  suggested  by  linear  finite  elements  (with  piecewise  linear  approximation 
of  the  source  distribution,  see,  e.g.,  [11,  pp.  45-46]) 

D.Ai  +  e  =  0  (6) 
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This  scheme  has  the  same  asymptotic  behavior  cis  the  pointwise  representa¬ 
tion,  but  improvement  in  performance  on  coarse  grids  is  evident  (see  Sec.  .3). 
For  variable  coefficients  this  becomes 


1  I  {k^<^)j+\ ^{k^9)j  +  (k^<l>)j-^  ,  /;4-i  +  4/j  + /j-i 
r0j  +  - - -  -I- - - - =  0  (7) 


where  ik^(t))j  = 


:i 
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On  a  non-uniform  grid  the  appropriate  weigliting  is 


3 


/ -  4>j  _  4>}  -  <f>}-i  \  j  d  h  ^ 

V  )/ 


0  (8) 


Superior  performance  on  non-uniform  grids  (see  See.  5)  is  attained  with  no 
increase  in  the  number  of  points  in  the  stencil.  For  a  discontinuity  in  physical 
properties  at  point  j  the  stetkcil  becomes 


I 

3 


( 


(l>]+\  -  4>i  <i>j-  4>3-\  \  //»■*■  +  h' 


h+ 


h~ 


+ 


0 


(9) 


Again,  this  may  also  be  considered  as  a  pieccwise-constant  approximation  of 
the  case  of  variable  coefficients. 

Performance  of  finite  difference  schemes  for  the  Helmholtz  equation  may 
be  enhanced  by  basing  the  stencils  on  more  general  defmitiojis  of  the  deriva¬ 
tive 


^  <A(x  4-  h)  - 
dx  0  li{kh)  h 


(10) 


where,  for  consistency 

lim  /?  =  1 

fch-O 


(11) 


This  definition  depends  on  fc/i,  an  indication  of  wave  resolution  by  the  grid. 
F'or  the  Laplacian  this  reduces  to  the  standard  definition  for  grids  of  any  size. 

This  generalization  of  the  derivative  definition  may  be  applied  for  either 
the  first  or  .second  derivatives,  or  to  both.  On  uniform  grids  all  are  equivalent. 
Since  the  parameter  dei)ends  on  the  grid  size  it  is  applicabh'  to  non-utiifonii 
grids  as  well.  Superior  performance  on  non-uniform  grids  is  obtained  by 


applying  this  concept  to  the  second  derivative  alone  (see  Sec.  5).  F'or  the 
uniform  case  this  reduces  to 


4>}+\  -  ^2  <^j+i  +  +  4>}-\  /j-n  +  ^fj  +  /j-i  _  y  ^  J2J 

/?/i*  '6  6 


The  resolution-dependent  parameter  0  is  defined  to  improve  method  per¬ 
formance.  For  example,  the  parameter  may  be  defined  to  eliminate  numerical 
dispersion 


6  1  —  cos{kh) 

{khy  2  cos{kfi) 


(1.3) 


so  that  in  simplified  settings  the  phase  is  exact  (EP),  resulting  in  no  trunca¬ 
tion  error  under  some  circumstances.  This  definition  satisfies  the  consistency 
requirement  (11).  In  such  cases  the  representation  of  source  terms,  which  is 
exact  for  piecewise  linear  variation,  is  no  longer  sufficiently  accurate.  A 
modification  of  the  representation  of  source  terms  that  does  not  degrade  the 
higher-order  accuracy  of  such  schemes,  similar  to  that  employed  by  HODIE 
methods  [2],  is 


~  -f  4>i-\  -I-  ^,-1  ^  /;  +  l/2  +  /j  4-  /j-1/2 


=  0  (14) 


(suggested  by  lineal  finite  elements  with  piecewise  quadratic  approximation 
of  the  source  distribution)  where  fj+\/2  is  the  load  term  evaluated  at  the 
midpoint  For  a  piecewise  linear  source  distribution  this  is  identical  to  (12). 
One  possibility  of  the  parameter 


/f  = 


12 

12  -  (kh)^ 


(15) 


yields  high-order  representation 

+  e  ^  '1*'  =  0  (16) 

(ho).  This  stencil  (without  the  modification  in  the  representation  of  source 
terms)  may  also  be  derived  by  employing  Fade  approximation 


1  -f  h^/\2Dr 


■4>,  -f  k^(f)j  +  fj  —  0 


(17) 


(see,  e.g.,  (12,  p.  538]). 

This  concept,  in  its  original  form,  which  may  be  viewed  as  an  average  of 
the  pointwise  and  weighted-average  representations  of  the  undifferentiated 
term  [13],  provides  high-order  performance  on  uniform  grids,  but  severely 
degrades  in  the  non-uniform  case.  However,  an  appropriate  generalization  to 
non-uniform  grids,  based  on  the  concept  of  generalizing  the  derivative  defi¬ 
nition,  leads  to  improved  performance  in  the  general  Ctise  as  well  (see  Secs.  5 
and  6).  Allowing  discontinuities  in  physical  coefficients  and  accounting  for 
non-uniform  grids  the  proposed  scheme  is 


( 


-  <f>j  4>j-\ N  -t-  /?■  h~ 


h+ 


h- 


-f 


1  /  {k-n-h-  \ 


K; 


/3+/i+ 


J+h+  +  0-h 
where  =  f3{k^h^). 


+  fj)  + 


0-h- 


0+h-^  +  0-h 


rTra  +  2/, 


i-1/2)] 


(18) 

0 


2.2  Two  dimensions 

Consider  a  two-dimensional  uniform  grid  of  size  h  with  points  at  ar,  =  ih 
and  j/j  =  jh.  For  simplicity  we  consider  the  homogeneous  case.  A  typical 
starting  point  is  the  five-point  representation 


"b  -f-  h 

^•.>+1  ~  2(^1, j  -t-  -2 

.  2  "  +  9i,j 


0  (19) 


which  is  the  two-dimensional  analog  of  (3).  Non-uniform  grids  and  material 
discontinuities  may  be  accounted  for  by  generalizations  of  (4)  and  (5). 

The  two-dimensional  counterpart  of  the  idea  that  leads  to  (6),  obtained 
by  employing  bilinear  finite  elements,  is  a  nine-point  representation 


+  2<pi^\,j  -t-  <^,4.i.j  +  i  —  -f  -b  2(j>i^\  j  -b 

+  2<^,,;+i  -b  —  80, .j  -b  +  20,.j_i  -I- 


-b 


6 


6 


(*^i+i  j+i  +  j  +  4^,-j+i  +  +  16^, ■j+ 


(20) 

=  0 


leading  to  a  significant  reduction  in  spurious  anisotropy  (see  Sec.  4).  HODIE 
methods  [2]  also  employ  nine-point  stencils  in  two  dimensions.  The  band¬ 
width  of  the  ensuing  linear  algebra  problem  is  typically  slightly  larger  but 
the  difference  in  the  cost  of  computation  is  insignificant. 

Performance  can  again  be  improved  by  substituting  fik  for  /i  as  in  the  one¬ 
dimensional  case  (14),  based  on  the  same  definitions  (13)  and  (15),  although 
the  methods  are  higher  order  only  for  propagatioti  along  grid  lines. 

In  order  to  maintain  higher-order  performance  on  uniform  grids  in  two 
dimensions  in  all  directions  of  propagation,  the  Fade  approximation  concept 
is  employed.  The  two-dimensional  counterpart  of  (17)  is 


1  + 


D. 


vv 


+  ~  D 
^  12 


(21) 


This  may  be  generalized  to 

4-  —  Dxx4>i,i  +  ^1  +  —  Dyycf) 
t’  (l  +  ^(D„  +  0„)  +  Tj^  D„D,\ 


+ 


=  0 


(22) 


where  7  is  selected  to  further  improve  properties  in  directions  other  than 
along  grid  lines,  without  effecting  dispersion  along  grid  lines  and  without 
degrading  higher-order  behavior  in  all  directions. 

The  standard  Fade  approximation  is  obtained  by  selecting  7  =  1  which 
yields  the  scheme 


( 

(  Dxx  +  +  —  Dzx  Dyy  I  4>i,}  + 

lO0i.j-l-l  +  +  100(^1, j-|- 

+  10<^,-,j_i  -f  10<^,_i,j  4 


(23) 

0 
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where 


yDrr  +  Dyy  +  —  DrxDyyj  0,,j  — 

4>i+\,j-\  +  S4>i+i.j  +  <Ai+i,_,+i  --  +  (ii,-\,j+i  I 

_ 

01-1, >+i  -f  —  2G(t>ij  -f  +  4>,+\j-\ 

1,2  ^  ^ 

Neglecting  higlier-order  terms  in  the  Fade  approximation  by  selecting  7  =  0 
leads  to  the  slightly  simplified  stencil  presented  in  [12,  p.  542) 


+  Dyy  +  —  DxxDyyj  0i,j  + 

+  4>ij+i  +  80,  j  +  0,,j-i  +  01-1, j)  =  0  (25) 

The  computational  cost  is  essentially  unaflfected  since  the  bandwidth  of  the 
algebraic  equations  is  identical.  Another  alternative  presented  in  [12.  p.  542] 


\ 

Dxx  +  Dyy  +  —  Dxx^yyj 


4"(0«+i.j+i  +  40,+i  j  +  40i,,+i  +  0i4.i,j_i  +  520i,_;+  (26) 

o 

0i-l.J  +  »  +  40i.J-t  +  40,-I,j  +  0,_l,;_t)  =  0 

obtained  by  selecting  7  =  2. 

Other  values  for  7  lead  to  other  alternatives.  In  Sec.  4  it  is  seen  that 
selecting  7  =  2/5,  which  leads  to  the  stencil 


Dxx  +  Dyy  +  —  Dxx  Dyy 

i  0 


0..r 


^(0i+i.j+i  +  2801+1,;  +  2S0,,j+i  +  0,+i,;_i  +  2440,,;+  (27) 

01-I.J+I  +280,,;_1  +280,-1,;  +0i-l,;-l)  =  0 

minimizes  dispersion  along  the  diagonals.  On  the  other  hand,  reducing  sen¬ 
sitivity  of  the  scheme  to  direction  of  propagation  is  attained  by  the  choice  of 
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7  =  14/5,  which  yields  the  stencil 

+  Dyy  4-  —  Dj-jrDyy^ 

P 

^(7<^,+i,j+i  +  16<^i+i,j  4-  16<^,,j4.|  +7(j>,+],j-\  +  261^4/, ,j+  (28) 

j+i  +  4-  lh0i_ij  4-  74>t-\.j~\)  =  U 

All  these  alternatives  reduce  to  (HO)  in  one  dimension.  Thus  the  dispersion 
analysis  for  (ho)  in  Sec.  3  describes  the  dispersion  of  all  alternatives  along 
grid  lines.  In  Sec.  4  the  performance  of  various  alternatives  in  other  directions 
is  compared. 


3  Spurious  Dispersion 

A  homogeneous,  isotropic  continuum  is  non-dispersive.  Waves  travel  at  a 
phase  velocity 

(29) 


UJ 

Cp  -  Co 

and  energy  propagates  at  the  group  velocity 

duj 


(30) 


and  so  both  are  identical. 

For  the  discrete  representation  this  is  usually  no  longer  the  case.  Each 
stencil  can  be  characterized  by  an  approximate  wav«*  number  =  k^{kh), 
which  depends  on  wave  resolution  and  thus  accounts  for  nnnn'rical  dispersion. 
The  pha.se  velocity  in  the  discrete  case  i.s 


u; 


k 


S  f,h  -  f,u 


(31) 


and  the  numerical  group  velocity  is 


= 


dijj 

()uj  Ok 

\  ()k 


C;i 


(32) 
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On  a  uniform  grid  in  one  dimension  a  numerical  solution  in  the  form  of 
a  plane  wave  is 

<l>j  =  exp{ik'^hy  (33) 

PT  For  point  j  the  pointwise  representation  (3)  of  the  plane  wave  solution 
yields 

0  =  e:xp{ik^h)  —  (2  —  {kh)^^  +  l/exp(ii^^) 

=  2cos{k^k)~  {2-ikhy)  (34) 

leading  to  the  dispersion  relation 

k^h  =  arccos  —  {khy/2^  (35) 

In  one  dimension  the  number  of  grid  points  in  a  wave  is 

G  =  27r /{kh)  (36) 

The  discrete  solution  represents  propagation  in  the  range  kh  <  2  which 
corresponds  to  a  limit  of  approximately  three  grid  points  per  wave¬ 
length.  Within  this  range  the  numerical  phase  velocity  is  thus 

Cp/co  =  kk/ axccos  —  {khy/2^  (37) 

and  the  numerical  group  velocity  is 

cj/cft  =  yJl-{khyiA  (38) 

Both  are  slower  than  the  speed  of  sound  in  the  material  cq. 


WA  Similarly,  for  the  weighted-average  representation  (6)  the  dispersion 
relation  is 


k^'h  —  arccos 


(\-{khy/^\ 

\l+{kh)y6) 


(39) 


representing  propagation  in  the  range  kh  <  \/T2,  a  limit  of  approxi¬ 
mately  two  grid  points  per  wavelength.  Within  this  range  the  numerical 
phase  velocity  is  obtained  directly  from  the  dispersion  relation  and  the 
numerical  group  velocity  is 


cj/co  =  y/l~{kh)y\2  (l  +  {khy/6)  (40) 

Both  are  faster  than  the  speed  of  sound  in  the  material. 
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EP  The  resolution-dependent  parameter  3  may  be  defined  so  that  discrete 
representations  are  non-dispersive  ( 13)  as  is  the  case  for  the  continuum. 
In  one  dimension  this  formulation  has  zero  local  truncation  error  for 
the  homogeneous,  constant  coefficient  case  on  uniform  grids,  and  the 
phase  and  group  velocities  are  exact.  Careful  generalization  leads  to 
improved  performance  on  general  configurations. 


HO  The  higher-order  representation  (15)  is  an  approximation  of  the  exact 
phcise  definition  (13)  on  uniform  grids.  The  resulting  higher-order  dis¬ 
persion  relation 


k^h  = 


arccos 


f\-b{khYl\2\ 

V  l+(a)Vl‘2  } 


(41) 


is  a  ( 1 , 1 )  Fade  approximation,  representing  propagation  in  the  range 
kh  <  \/6,  a  limit  of  approximately  2  1/2  grid  points  per  wavelength. 
Within  this  range  the  numerical  phase  velocity  is  again  obtained  di¬ 
rectly  from  the  dispersion  relation  and  the  numerical  group  velocity 

is  _ 

c)lcc  =  sJT-  {khfl^  (l  4-  {khfin)  (42) 

Both  are  slower  than  the  speed  of  sound  in  the  material.  The  power 
series  expansion  of  the  dispersion  relation 


k^'h 


(H)’ 


480  12096 


>  kh 


(43) 


demonstrates  the  higher-order  nature  of  this  representation. 


Dispersion  of  the  various  formulations  is  plotted  in  Fig.  1 .  Note  that  the 
region  of  primary  interest  is  6’  >  4,  a  resolution  of  at  least  four  grid  points 
per  wavelength.  Within  this  region  the  errors  in  the  pointwise  and  weighted- 
average  representations  are  similar,  and  the  asymptotic  behavior  is  the  same 
(see  Sec.  6).  However,  even  approaching  the  limit  of  resolution,  and  certainly 
beyond  it,  the  weighted-average  performance  is  superior.  For  example,  at 
the  limit  of  resolution  (G  =  4)  there  is  a  38%  error  in  the  group  velocity  for 
the  pointwise  representation,  whereas  in  the  weighted-average  rejjresentation 
the  error  is  only  26%.  The  higher-order  method  offers  significantly  superior 
representation,  an  error  of  only  7%  in  group  velocity  at  G  =  4,  and  the 
exact-phase  method  provides  dispersion-free  solutions. 


Group  velocity,  c  /c  Phase  velocity. 


4  Spurious  Anisotropy 

On  the  uniform  grid  in  two  dimensions  a  numerical  solution  in  the  form  of  a 
plane  wave  oriented  at  angle  6  to  the  grid  lines  is 

h co^  Oy  exp{ih*^ h  sm  Oy  (44) 

PT  For  point  i,  j  the  pointwise  repn^entation  (19)  of  the  plane  wave  solution 
yields 

0  =  exp(t^‘^/i  sin  9)  +  exp{ik^h  cos  0)  —  ^4  —  (kh)^^  + 

I  /  exp{ik^k  cos  0)  +  1/  exp{tk^h  sin  0) 

=  2  (cos{k^h  cos  6)  +  cos(A:^*/i  sin  0))  —  ^4  —  {khy^  (45) 

leading  to  various  dispersion  relations,  depending  on  the  angle  of  orien¬ 
tation  6.  When  the  wave  is  aligned  with  the  grid  (e.g.,  0  =  0)  this  leads 
to  the  one-dimensional  dispersion  relation  (35).  The  other  extreme  case 
occurs  when  the  wave  is  oriented  in  the  direction  of  cell  diagonals  (e.g., 
9  =  7r/4) 

k^h  =  \/2arccos  ^1  —  {kh^/Aj  (46) 

The  discrete  solution  represents  propagation  in  the  range  kh  <  \/S. 
Within  this  range  the  numerical  phase  velocity  is  again  obtained  di¬ 
rectly  from  the  dispersion  relation  and  the  numerical  group  velocity 
is 

cj/co  =  \/l  -  {khy/S  (47) 

Both  are  slower  than  the  speed  of  sound  in  the  material.  It  is  interesting 
to  note  that  the  pointwise  representation  is  more  dispersive  when  waves 
are  oriented  along  the  grid. 

WA  Similarly,  for  the  weighted-average  representation  (21)  dispersion  rela¬ 
tions  are  obtained  from 

4  ^6  —  {khy^  —  ^12  -I-  {khy^  cos(A:^/j  cos  0)  cos(fc^/i  sin  0)—  (48) 

2  ^3  -I-  {kh)^'\  (coii(k''h  cos  0)  -h  cos{k^h  sin  0))  = 
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t) 


When  the  wave  is  aligned  with  the  grid  this  leads  to  the  known  one¬ 
dimensional  dispersion  relation  (39)  and  for  waves  parallel  to  cell  diag¬ 
onals 


k''h  =  \/2arccos 


( 1  -  mya  \ 

[l-h(kWl2j 


(49) 


representing  propagation  in  the  range  kk  <  y/24.  Within  this  range 
the  numerical  ph£ise  velocity  is  obtained  directly  from  the  dispersion 
relation  and  the  numerical  group  velocity  is 


cj/co  =  sjx  -  (fc/i)V24  (l  -f  [khfin]  (50) 

Both  are  faster  than  the  speed  of  sound  in  the  material.  The  weighted- 
avtrage  representation  is  also  more  dispersive  for  waves  that  are  ori¬ 
ented  with  the  grid. 


HO  Dispersion  relations  for  the  higher-order  representation  with  ft  defined 
in  (15)  at  various  angles  of  orientation  may  be  found  in  similar  fashion. 
For  waves  aligned  with  cell  diagonals  the  relation  is  identical  to  that  of 
the  pointwise  representation  at  this  angle  (46),  and  the  same  holds  for 
the  wave  velocities.  This  representation  is  thus  higher  order  only  for 
waves  oriented  along  grid  lines. 

Representations  that  ere  truly  higher  order  in  all  directions  of  propa¬ 
gation  are  based  on  (22).  For  waves  aligned  with  the  grid  this  leads  to 
the  higher-order  dispersion  relation  (41)  and  along  cell  diagonals  the 
relation  is 


k^h  =  v^arccos  f^V  ^  "  7)(^-MVl44  -  (4  -K6  -  i){kh)^l\2) 


2-|-7(fc/i)Vl2 

By  examining  the  power  scries  expansion  of  this  relation 


(51) 


k^h  Ks  kh  -X-  (57  -  2)^ — -  -f  — — 
'  '  ’  5760  96768 


(52) 


it  is  clear  that  the  value  of  7  =  2/5  minimizes  dispersion  in  the  direction 
of  cell  diagonals.  On  the  other  hand,  7  =  14/5  minimizes  the  difference 
between  the  dispersion  along  grid  lines  and  in  the  direction  of  cell 
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diagonals,  as  seen  by  comparison  to  the  power  series  expansion  of  the 
higher-order  dispersion  relation  along  grid  litics  (4il).  Dispersion  in  the 
direction  of  cell  diagonals  for  various  values  of  7  is  plotted  in  Fig.  2.  As 
expected,  the  stencil  with  7  =  2/5  is  essentially  tion-dispersive  in  the 
region  of  primary  interest  with  a  resolution  of  at  least  four  grid  points 
per  wavelength. 

The  ratio  between  numerical  iispeision  along  grid  lines  and  along  cell 
diagonals  is  shown  in  Fig.  3.  By  design,  the  stencil  with  7  =  14/5  is  the 
least  anisotropic  in  the  range  of  at  least  four  grid  points  per  wavelength. 
This  is  corroborated  by  the  polar  plots  in  Fig.  4,  showing  the  variation 
in  phase  velocity  with  angle  of  orientation  for  various  resolutions.  Note 
that  the  figure  shows  to  accentuate  deviations  from  the  exact 

value  of  unity.  As  expected,  all  the  schemes  perform  identically  along 
grid  lines,  but  behavior  in  other  directions  is  determined  by  the  choice  of 
7.  Differences  among  the  various  cases  become  more  pronounced  with 
reduced  resolution,  but  in  general  are  not  extreme.  The  two  schemes 
that  stand  out  are  indeed  7  =  2/5  which  minimizes  dispersion  along 
diagonals,  and  hence  overall,  and  7  =  14/5  which  redtices  anisotropy. 

EP  The  case  that  is  non-dispersive  in  one  dimension  (13)  may  be  treated 
similarly.  In  this  case  there  is  no  dispersion  for  waves  aligned  with 
the  grid  and  the  dispersion  relation  for  waves  in  the  direction  of  cell 
diagonals  is 

tfxL  (  \+2cos{kh)\ 

k  h  =  v2arccos  (  2  — - -  ---  (5.1) 

\  5  +  cos(il-/i)  /  ^  ’ 

The  numerical  pheise  velocity  is  again  obtained  direcily  from  the  dis¬ 
persion  relation  and  the  numerical  group  velocity  is 


7-t-5cos(il-/i) 
1  -f  cos(kh) 


5  -f-  cos(A:/t)'\ 

,  ) 


This  representation  is  obviously  less  dispersive  for  waves  that  are  ori¬ 
ented  with  the  grid. 

Dispersion  in  the  direction  of  cell  diagonals  of  the  various  formulations 
is  plotted  in  Fig.  5.  Recall  that  the  region  of  primary  interest  is  G  >  4,  a 
resolution  of  at  least  four  grid  points  per  wavekmglh.  Dispersion  properties  of 


Phase  velocity. 


Group  velocity 


Figure  4:  Polar  plots  of  anisotropy  in  (ep/ro)'*  of  higher-order  tliscrete  for¬ 
mulations  based  on  generalized  Fade  approximation. 
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each  scheme  at  arbitrary  orientations  are  bounded  on  one  hand  by  dispersion 
along  grid  lines  shown  in  Fig.  1,  and  on  the  other  hand  by  dispersion  along 
cell  diagonals  shown  in  Fig.  5.  Performance  of  the  pointwisc  and  weighted- 
average  schemes  improves  the  farther  the  orientation  of  propagation  is  from 
the  direction  of  grid  lines.  The  same  holds  for  the  higher-order  .scheiiK's  of 
interest  with  7  =  2/.')  and  7  =  14/5.  For  the  (EP)  method  the  opposite 
occurs,  so  that  performance  of  this  scheme  is  vastly  superior  along  grid  lines. 
This  scheme  is  higher  order  only  along  grid  lines,  but  it  still  maintains  a  high 
degree  of  phase  accuracy  in  all  orientations. 

The  resolution-dependent  parameter  (i  may  be  defined  so  that  the  nu¬ 
merical  representation  is  non-dispersive  for  waves  at  any  given  angle  of  ori¬ 
entation.  For  example 

„  o  1  -  cos  (•/2kh/2) 

(UP  2  +  cos  (,/2W./2)  '■ 

eliminates  dispersion  of  waves  along  cell  diagonals.  Similar  performance  was 
attained  in  the  context  of  finite  element  methods  [14].  In  general,  how¬ 
ever,  the  direction  of  wave  propagation  is  not  known  in  advance  and  there 
is  a  concern  that  defining  /?  for  any  orientation  other  than  along  grid  lines 
may  degrade  performance  on  non-uniform  grids,  as  discussed  in  the  follow¬ 
ing  section.  Furthermore,  grids  should  be  aligned  with  dominant  directions 
of  propagation  to  the  extent  possible.  For  these  reasons  it  is  preferred  to 
maintain  dispersion-free  discrete  solutions  along  grid  lines. 

Numerical  dispersion  is  thus  sensitive  to  the  orientation  of  wave  propaga¬ 
tion.  The  two  extreme  cases  are  along  grid  lines  shown  in  Fig.  1,  and  along 
cell  diagonals  shown  in  Fig.  5.  The  largest  change  in  dispersion  properties 
possible  is  thus  the  ratio  between  the  two,  shown  in  Fig.  6.  Recall  that  the 
region  of  primary  interest  is  G  >  4,  a  resolution  of  at  least  four  grid  points  per 
wavelength.  For  highly  resolved  phenomena  the  performance  of  all  schemes 
is  similar  and  quite  good.  As  wave  resolution  is  reduced  only  the  higher- 
order  schemes  (with  values  of  7  shown)  retain  a  low  level  of  anisotropy.  Of 
the  other  schemes,  approaching  the  limit  of  resolution  and  certainly  beyond 
it,  the  pointwise  method  is  clearly  more  sensitive  to  direction  of  propaga¬ 
tion.  For  example,  at  the  limit  of  resolution  (6'  =  4),  there  is  8%  anisotropy 
in  the  pointwise  representation  of  phase  velocity,  whereas  the  anisotropy  of 
other  methods  is  at  most  about  half  that  value.  This  becomes  even  more 


.  5;  Phase  ami  group 
cell  tliagona's- 


velocities 


of  iwo-d'nmMisio 


iial  dibcrele 


pronounced  in  group  velocity. 

Figure  7  shows  the  variation  in  phase  velocity  with  angle  of  orientation 
of  different  schemes  at  various  resolutions.  For  presentation  purposes  the 
figure  shows  (Cp/co)'*.  Note  that  the  plot  for  (PT)  does  not  include  the  Ccise 
of  G  =  3  since  this  scheme  no  longer  represents  propagation  at  this  low 
resolution,  which,  in  any  event,  is  outside  the  region  of  primary  interest  of 
G  >  4.  It  is  clear  from  these  plots  that  the  numerical  phase  velocity  is  less 
than  the  speed  of  sound  in  the  material  in  all  cases  shown  except  for  (WA). 
C’lose  examination  of  Fig.  4  indicates  that  this  is  true  of  higher-order  methods 
only  with  7  >  2/5.  With  the  exception  of  (EP),  all  the  schemes  considered 
exhibit  superior  dispersion  behavior  along  cell  diagonals.  This  would  not  hold 
for  higher-order  methods  with  7  >  14/5,  but  there  is  no  apparent  motivation 
to  pursue  such  methods  in  the  first  place.  As  mentioned,  employing  (55) 
eliminates  dispersion  along  cell  diagonals,  leading  to  a  version  of  (ep)  with 
superior  dispersion  behavior  along  diagonals  that  is  similar  to  other  schemes 
in  this  regard. 

Overall,  high  wave  resolution  or  higher-order  methods  are  required  if 
anisotropy  is  a  concern.  Of  the  methods  that  are  not  high  order,  on  grids  with 
lower  resolution,  the  weighted  average  representation  and  its  enhancements 
are  much  less  anisotropic  than  the  standard  pointwise  representation. 

5  Spurious  Reflection  and  Transmission 

Reflected  and  transmitted  waves  are  generated  by  incident  waves  at  disconti¬ 
nuities  in  physical  properties.  Numerical  di.spcrsion  of  discrete  formulations 
gives  rise  to  incorrect  representation  of  these  phenomena  at  transitions  in 
wave  resolution. 

5.1  Grid  transition 

In  a  homogeneous  material  no  reflection  should  occur.  However,  changes  in 
grid  size  alter  wave  resolution  giving  rise  to  spurious  reflection  and  transmis¬ 
sion  due  to  numerical  dispersion,  phenomena  that  may  be  characterized  in 
a  manner  similar  to  that  of  waves  at  discontinuities  in  physical  coefficients. 
The.se  phenomena  are  well  known  [15],  have  been  carefully  analyzed  [16]  and 
numerically  demonstrated  [17]. 
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Figure  6:  The  ratio  of  numerical  dispersion  along  grid  lines  and  along  cell 
diagonals  for  two-dimensional  discrete  formulations. 


(’oiisider  one-dinionsional  configuration  discretized  by  a  unifonn  grid 
of  size  h~  to  .ie  left  of  the  origin  and  a  uniform  grid  of  size  /i"*"  to  its  right. 

Due  to  the  transition  in  grid  size  at  the  origin,  an  incident  plane  wave  of  unit 

niagiiitude  traveling  in  the  positive  direction  with  discrete  values  given  in 
(33)  for  j  <  0  generates  spurious  transmission  such  that  </»o  ^  1  and  spurious 
reflection  of  magnitude  <i>o  —  \,  The  numerical  solution  is  thus 

where  the  dispersion  error  is  represented  by  the  numerical  wave  numbers 
it'**  =  k^{kh^),  and  the  transmission  error  is  represented  by  <f>o.  In  particular 

<A,  =  <l>oexp{ik''*  h*)  (57) 

4>-i  =  ^oexp(zA:*  h~)  —  2is\n{k^  h~)  (58) 


PT  The  pointwise  representation  (4)  of  this  solution  at  the  origin  yields 

2  /■<>,-  (1  -  (H+)V2)«o  .  ^-1  -  (1  -  (H-)V2W„'\ 

“  "  F+T^v  ie  *  F  ) 

2  /  4>oc^p{ik^*  h'*')  —  cos(fc^'’^/i‘^)(^o 

“  A+T/i-  V  “*■ 

<t>oexp{ik^~ h~)  —  2i  s\n{k''~ h~ )  ~  cos{k''~  h~)<f>o^ 


where  the  second  line  was  obtained  by  the  dispersion  relation  for  the 
pointwise  representation  (35).  Thus 


2yJ\-{kh-YI^ 

^\-(kh-Yi\^^\-[kh^YlA 


(60) 


which  is  valid  in  the  range  of  resolution  in  which  the  pointwise  formu¬ 
lation  represents  propagation  (along  gri'l  lines). 
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WA  Similarly,  the  weighted-average  representation  (8)  of  the  solution  (56) 
at  the  origin  leads  to  spurious  transmission 


<i>Q  = 


2^\-{kh-YI\2 


^1  -(M-)Vl2  +  0  -(A:/i+)Vl2 


(61) 


which  is  valid  in  the  range  of  resolution  in  which  the  weighted-average 
formulation  represents  propagation  (along  grid  lines). 


HO  Transmission  for  the  higher-order  representation  with  the  parameter 
defined  in  (15)  is 


4>q  = 


{\-{kh-YI\2) 


^\-ikh~)ye  ^i-{kh+Y/6 

(1  -  {kh-y/\2)  (1  -  {kh*Y/\2) 


(62) 


which  is  valid  in  the  range  of  resolution  in  which  the  higher-order  for¬ 
mulation  represents  propagation  (along  grid  lines). 

EP  The  case  that  is  non-dispersive  in  one  dimension  (13)  may  be  treated 
similarly.  In  this  case  the  transmission  is 


^  sin(A;fe  )/h 

_ 2  -h  cos{kh~) _ 

s\n{kh~)/h~  s\n{kh'^)/h'^ 

2  -|-  cos(kh~  )  "^^  2  -I-  cos{kh'^) 


(63) 


which  is  valid  in  the  range  of  resolution  in  which  the  higher-order  for¬ 
mulation  represents  propagation  (along  grid  lines). 


Spurious  transmission  of  the  various  formulations  at  different  wave  resolu¬ 
tions  is  plotted  in  Fig.  8.  In  general  the  sensitivity  to  transition  in  grid  size  is 
higher  for  coarser  grids.  The  weighted-average  representation  is  significantly 
siiperior  to  the  pointwise  scheme  on  non-uniform  grids.  The  higher-order 
and  exact-phase  formulations  offer  further  improvement. 
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Transmission  error  (G 


5.2  Interface  of  physical  properties 


Discontinuities  in  physical  properties  give  rise  to  wave  reflection  and  trans¬ 
mission.  The  relative  amplitudes  of  the  reflected  and  transmitted  waves 
depend  on  the  ratio  of  wave  numbers,  which  defines  the  character  of  the 
discontinuity.  The  numerical  representation  of  these  phenomena  by  finite 
element  methods  was  studied  in  [18). 

Consider  a  generalization  of  the  previous  configuration  in  which  a  discon¬ 
tinuity  in  material  properties  as  well  as  a  jump  in  grid  size  may  occur  at  the 
origin,  so  that  k~  is  the  wave  number  to  the  left  of  the  origin  and  — to 
its  right.  An  incident  plane  wave  of  unit  magnitude  traveling  in  the  positive 
direction  exp(tA;~x)  for  x  <  0  generates  reflected  and  transmitted  waves,  so 


that 


(  exp(tfc“x)  -f  (^(0)  —  1)/  exp(tA:‘'‘x),  x  <  0 
^  <^(0)  exp(tfc‘*'x),  X  >  0 


(64) 


where 


m  = 


2  k- 

k-  f  k+ 


(65) 


The  discrete  solution  is  again  (56)  where  the  numerical  wave  numbers  are 
=  k^^ikH^). 


PT  The  pointwise  representation  (5)  of  the  solution  at  the  origin  yields 


0  = 


h+  +  h 
2i 


-( 


A+)V2)^o  ,  ^-^-{\-{k-h-f|2)<i>, 


h+ 


-f 


s\n{k^* h'*')  ^  s\n{k^~h  ) 


h+ 


h- 


(f>o  —  2- 


h- 

s\n{k^~  h-y 


h- 


(66) 


where,  again,  the  dispersion  relation  for  the  pointwise  representation 
(35)  is  employed.  Thus 


2k-,Jl-{k-h-y/4 

k-^\  -{k-h-)y4  +  k+^\  -  [k+h+yi4 


(67) 
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Physical  transmission  depends  on  the  ratio  of  the  wave  numbers.  Nu¬ 
merical  solutions  depend  on  this  parameter  and  on  the  ratio  of  re.solutions. 
To  find  out  which  of  the  two  parameters  significantly  effects  the  numerical 
error  in  transmission  consider  the  transmission  error  as  a  function  of  ratio  of 
resolutions  for  6  grid  points  per  wavelength  to  the  left  of  the  origin.  This  is 
plotted  for  a  ratio  of  the  wave  numbers  equal  to  unity  in  Fig.  8  (top).  Increeis- 
ing  the  ratio  by  one  order  of  magnitude  and  by  two  yields  the  behavior  shown 
in  Fig.  9.  The  difference  between  these  plots  is  not  significant  indicating  that 
the  error  depends  primarily  on  the  ratio  of  resolutions.  All  the  representa¬ 
tions  have  the  property  that  (fiQ  =  ^(0)  if  k'^h'^  =  k~h~.  Again,  superior 
performance  of  the  weighted-average  representation  and  its  enhancements  is 
evident. 
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Figure  9:  Transmission  error  (4>o  —  ^(0))/<^(0)  at  wave  number  ratios  of  10 
(top)  and  100  for  6  grid  pomts  per  wavelength  to  the  left  of  the  origin. 
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6  Local  Truncation  Error  Analysis 

The  local  truncation  error  is  the  residual  left  by  substituting  the  exact  solu¬ 
tion  in  the  discrete  representation.  In  the  following,  sufficient  differentiability 
is  8issumed  for  all  functions  involved. 


6.1  Uniform  grids 

Consider  the  one-dimensional  constant-coefficient  case  on  a  uniform  grid.  For 
the  pointwise  representation  (3) 

^  ^  ^  ^ 

=  rM  +  ^r(rj) + (*■") + ^"‘(^■))  +  *V(x>) + /(x,) 

=  ^r(xj)  +  o(h‘)  (71) 

where  primes  and  superior  Roman  numerals  indicate  differentiation  by  the 
argument.  The  second  line  is  obtained  by  Taylor’s  formula,  where  > 
x"  >  Xj  and  Xj  >x'*‘>  Xj+i ,  and  the  third  line,  which  follows  from  the  fact 
that  0  satisfies  the  Helmholtz  equation,  indicates  consistency.  The  pointwise 
scheme  is  thus  second-order  accurate. 

The  weighted-average  representation  (6)  is  similar  on  uniform  grids 

^(Xj+i)-2(^(xj)-f  ,  ,3  (^(x>+i) +  4<^(xj)-|- ^(Xj_,)  , 

T  ~  - — - 1-  fc  - ^ -  -f 

/(X,>l)-m(X;)-f  /(x,-l) 

=  (l  +  4  + 

J[Xj)  -f  - ^ - 

=  +  o(h*)  (72) 

The  weighted-average  scheme  is  also  second-order  accurate.  This  order  of 
accuracy  is  retained  in  the  case  of  variable  coefficients  (7). 
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For  the  improved  representations  (14) 


^  6 

f{^j+\/2)  +  f{Xj)  +  f{Xj-l/2) 

/(^j-t-i/i)  -  2/(xj)  +  /(jj-1/2) 


Employing  the  definition  of  p  that  leads  to  the  high-order  representation  (15) 
yields 

r  =  -  rMW  +  0(h‘)  (74) 

justifying  its  name  as  a  higher-order  scheme.  Note  that  if  the  source  terms 
were  not  represented  appropriately  there  would  be  second-order  terms  in  the 
truncation  error.  The  scheme  that  is  dispersion-free  in  one  dimension  (13) 
has  a  truncation  error 

r  =  -  r(i;)/4) + o(a»)  (75) 

If  the  fourth  derivative  of  the  source  vanishes  the  method  becomes  six-order 
accurate.  Furthermore,  the  truncation  error  is  zero  when  all  the  derivatives 
of  the  source  from  fourth  order  and  higlier  vanish. 


6.2  Non-uniform  grids 

In  analyzing  method  performance  on  non-uniform  grids  a  change  of  variables 
from  physical  space  to  computational  space  is  often  considered,  so  that  the 
grid  is  uniform  in  the  latter  [19],  The  order  of  accuracy  of  some  methods 
on  non-uniform  grids  may  drop  in  physical  space.  Nevertheless,  in  computa¬ 
tional  space  it  remains  unchanged  from  the  order  in  the  uniform  case.  To  a 
certain  extent  grid  stretching  should  reflect  variation  of  the  solution.  In  this 
case  accuracy  in  computational  space  is  representative  of  the  situation.  In 
practice,  however,  grid  variation  is  determined  by  geometric  considerations 
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as  well.  The  following  results  are  thus  presented  in  physical  space,  which 
describes  the  more  general  case. 

For  the  pointwise  representation  (4) 


T  = 


( -  <i>{xj)  /h-*-  H-  h- 


I 


h+ 


h- 


2 


+  +  fixj) 


-  <f>  (Xj)  + - - - (f)  {Xj)  + - ^  \^J>  + 


1 


60(/t+  +  +  fix,) 


-h- 


[h^Y  -h*h-  ^{h-y 
12 


r(x,)+ 


r{x,)  +  o{e) 


(76) 


The  pointwise  scheme  is  indeed  second-order  accurate  on  uniform  grids,  but 
may  drop  to  first  order  in  the  non-uniform  case. 

Whether  the  scheme  actually  drops  to  first  order  or  not  depends  on  the 
degree  of  grid  stretching.  If 


=0(/i»’+'),  p>0  (77) 


the  stretching  is  called  algebraic  [19].  With  algebraic  stretching  the  pointwise 
scheme  retains  second-order  accuracy.  Otherwise  the  accuracy  drops  to  first 
order. 

In  contrast,  for  the  weighted-average  representation  (8) 


r  = 


( 0(3^j+i)  -  ^(xj)  _  0(xj)-d>(x,-i)\  /h^  +  h 
\  h-  )/  2  ^ 

T  (/t+ +  +  <t>{xj-\)'^  + 


1  ( 

3  \h+  +  h- 


h- 


h+  +  h 


<^(Xj>i)  -  4>(Xj) 

6  J  h+ 

{kh-Y\  4>(x,.,)  -  <f>{x 


3(2/(i,)  +  /(i 


j-.))) 


r  ((l  +  I  -■ 


( 


1  + 


h- 


— ^  +  k'^(f>{Xj 


)  + 
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first-order  terms  cancel  out  even  on  non-uniform  grids  so  that  second-order 
accuracy  is  retained  in  any  case.  Again,  these  results  apply  to  the  ceise  of 
variable  coefficients  cis  well. 

For  the  improved  representations 

’■  =  [ — hf - r- — )/ — 2 —  + 

-  2  // 

0*k*  +  ^-k-  V\  6  j  k* 

+r/>-(/(x.)-/(x,.,„)))  (79) 

For  the  definition  of  fi  that  leads  to  the  high-order  representation  (15) 


r  = 


{h^-h-)[{h^Y^{h-f) 


90 


(fcV"'(x,)  -  r(x,)/4)  + 


{h^Y  -  [h^fh-  -H  {h^  ?{h-f  -  h*{h-f  +  (/i-)\,  2 


240 


0{h^) 


(fcV’'(x;)-r(.r,)/4)  + 

(80 


the  truncation  error  is  third-order  accurate  on  non-uniform  grids  (and,  of 
course,  fourth  order  in  the  uniform  case). 


7  Conclusions 


In  this  work  finite  difference  methods  for  solving  problems  of  time-harmonic 
acoustics  are  developed  and  analyzed.  The  well-known  pointwise  represen¬ 
tation,  eqs.  (3)  and  (19)  in  one  and  two  dimensions,  respectively,  is  second- 
order  accurate  on  uniform  grids.  However,  accuracy  may  drop  to  first  order 
in  the  non-uniform  case  (4)  unless  sufficiently  smooth  grid  stretching  (77) 
is  employed.  In  multi-dimensional  configurations  the  representation  actually 
improves  the  less  aligned  the  propagation  directions  are  with  respect  to  the 
grid. 

A  weighted-average  representation,  eqs.  (6)  and  (21)  in  one  and  two  di¬ 
mensions,  respectively,  has  the  same  asymptotic  behavior  on  uniform  grids, 
but  is  less  sensitive  to  low  wave  resolution  and,  more  importantly,  to  di¬ 
rection  of  propagation  and  transition  in  wave  resolution  (including  material 
interfaces).  Performance  in  multi-dimensional  configurations  again  improves 
for  propagation  directions  that  are  not  aligned  with  the  grid.  In  general, 
anisotropy  in  numerical  representation  is  reduced  with  increased  wave  reso¬ 
lution.  At  lower  resolution  the  weighted-average  representation  (21)  is  much 
less  anisotropic  than  the  standard  pointwise  representation  (19).  Second- 
order  accuracy  is  retained  on  any  non-uniform  grid  (8)  at  virtually  no  in¬ 
crease  in  computational  cost.  These  results  hold  for  variable  coefficients  as 
well. 

Superior  performance  is  attained  by  basing  the  schemes  on  a  generalized 
definition  of  the  derivative  (10)  which  incorporates  a  resolution-dependent 
parameter.  Improved  schemes  with  higher-order  accuracy  are  designed  by 
appropriate  definition  of  the  parameter  (15),  reducing  spurious  dispersion 
and  reflection.  Defining  the  parameter  for  schemes  which  are,  in  some  cases, 
dispersion-free  (13)  leads  to  the  same  asymptotic  behavior  with  improved 
coarse  grid  accuracy.  Source  terms  must  be  represented  accordingly  (14)  so 
as  not  to  degrade  the  higher-order  accuracy.  These  methods  are,  in  general, 
fourth-order  accurate  on  uniform  grids  and  third  order  in  the  non-uniform 
case.  The  performance  of  these  schemes  in  multi-dimensional  configurations 
is  superior  for  any  direction  of  propagation.  Their  performance  improves 
as  propagation  directions  become  aligned  with  the  grid.  In  principle,  grids 
should  thus  be  aligned  with  directions  of  propagation  to  the  extent  possible, 
further  enhancing  the  performance  of  these  methods. 

Schemes  that  exhibit  higher-order  behavior  on  uniform  grids  in  all  di- 
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rections  in  two-dimensional  configurations  are  derived  on  the  b^lsis  of  Fade 
approximation  and  its  generalization  (22).  The  dispersion  of  these  methods 
(^ls  well  as  their  spurious  reflection  and  transmission)  along  grid  lines  is  iden¬ 
tical  to  that  of  the  higher-order  method  based  on  the  generalized  definition 
of  the  derivative.  Dispersion  along  grid  diagonals  is  minimized  by  employing 
7  =  2/5  which  leads  to  (27).  These  methods  are  by  far  less  anisotropic  than 
all  other  schemes.  The  value  of  7  =  14/5  leads  to  the  stencil  with  the  lowest 
degree  of  anisotropy  (28). 

In  general,  wave  resolution  (kh)  should  be  kept  as  even  as  possible  through¬ 
out  the  grid  to  minimize  spurious  reflection  and  transmission.  Sensitivity  to 
these  phenomena  is  greater  on  relatively  coarse  grids. 
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